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Abstract 

A  second-order  accurate  method  for  solving  the  combined  potential  and  circuit 
equations  in  an  electrostatic  bounded  plasma  PIC  simulation  is  presented.  The 
boundary  conditions  include  surface  charge  on  the  electrodes,  which  are  connected 
to  a  series  RLC  circuit  with  driving  terms  V(t)  and  I(t).  The  solution  is  obtained  for 
planar,  cylindrical,  and  spherical  electrodes.  The  result  is  a  tridiagonal  matrix  which 
is  readily  solved  using  well-known  methods.  The  method  is  implemented  in  the  codes 
PDP1  (Plasma  Device  Planar  ID),  PDC1  (Cylindrical),  and  PDS1  (Spherical). 

I.  INTRODUCTION 

A  comprehensive  review  of  the  considerations  involved  in  bounded  plasma  particle  simu¬ 
lation  is  presented  by  W.  S.  Lawson  [1].  Lawson  presents  a  method  which  is  second-order  accurate 
when  A t2/LC  <t  1  and  R  AtiL  «  1,  and  is  stable  for  A t2ILC  <  2  and  RAt/L  <  2.  Here  we  improve 
on  the  accuracy,  stability,  and  simultaneity  of  the  solution  for  potentials  in  a  bounded  one¬ 
dimensional  plasma  with  external  circuit  and  driving  terms. 

In  [1]  and  [2]  the  boundary  conditions  are  decoupled  from  the  potential  equation.  A  first- 
order  circuit  solution  is  used  when  the  inductance  is  zero.  The  scheme  is  self-consistent  when  L 
is  non-zero  and  the  applied  (driving)  potential  is  small  compared  to  the  space-charge  potential 
across  the  system.  These  conditions  are  violated  for  a  large  class  of  problems,  including  capac- 
itively  coupled  RF  discharges  and  plasma  immersion  ion  implantation  materials  processing; 
therefore,  a  new  method  is  desired. 

Particle-in-Cell  (PIC)  methods  weight  particles  to  a  spatial  grid  using  a  particle  shape  factor 
to  obtain  charge  and  current  densities  on  the  grid  [2].  For  example,  the  code  PDP1  uses  the  linear 
weighting  scheme  shown  in  Figure  1.  The  field  and  circuit  solution  presented  here  is  independent 
of  the  weighting  scheme  used;  we  assume  that  the  charge  density  is  given  on  the  spatial  grid. 
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single  panicle  density 


Figure  1 .  PDP1  PIC  with  linear  weighting  to  the 
spatial  grid.  The  subscript  / is  the  particle  index,  j is 
the  grid  index.  For  particles  in  a  cell  adjacent  to  an 
electrode,  the  weighting  must  also  account  for  the 
half  width  of  the  cell. 


Figure  2.  Particle  weighting  to  a  radial  grid,  using 
area  (cylindrical  model)  or  volume  (spherical  model) 
ratio. 


Particles  of  finite  size,  cylindrical  shells  in  PDC1  and  spherical  shells  in  PDS1,  are  placed 
in  a  gridded  system,  and  weighted  to  the  grid  to  obtain  p(ry)  at  the  grid  points.  The  particles  are 
assumed  to  have  uniform  density,  so  the  area  of  rings  or  the  volume  of  shells  can  be  used  to  weight 
the  charges  to  the  grid  as  shown  in  Figure  2.  The  particle  of  finite  width  Ar  is  centered  at  r,.  The 
intersection  between  the  finite  particle  and  the  grid  cell  determines  the  fraction  of  the  particle 
charge  assigned  to  each  grid  node  rr  This  is  cloud-in-cell  weighting,  versus  the  particle-in-cell 
weighting  in  [2,  Figure  14-1  la].  The  fractions  of  the  charge  assigned  to  the  grids  are 
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In  the  spherical  model,  the  squared  terms  are  replaced  by  cubic  terms. 

The  charge  density  on  the  grid  is  used  to  solve  the  Poisson  equation,  V20  =  -p/e,  or  the 
equivalent  flux  conserving  potential  equation  obtained  from  Gauss’  Law.  Once  the  potentials  are 
known,  the  electric  field  on  the  grid  can  be  obtained  from  E  =  -VO.  The  force  on  the  particles  is 
obtained  by  interpolating  E  at  each  particle  position  using  some  weighting.  The  particle  velocity 
is  updated  using  the  Lorentz  force  equation. 


dv, 

dt 


—  [ECqO  +  v.xB]  . 
m, 


(3) 


The  particle  positions  are  updated  using 


d\ 

dt 


(4) 


In  cylindrical  coordinates,  Eq.  (4)  is 


dvr  „  mvl 

m  — —  =  Fr  + - 

dt  r 


d(rv6) 

m  — =  rf8, 


and 


(5) 


Here  r,  0,  and  z  are  the  cylindrical  coordinates  and  vr  =  dr/dt  and  v0  =  rdQ/dt.  For  the  cylindrical 
model,  B  =  Btz  so  Fr  -  q(Er  +  F&  =  -qvJB„  and  F,  =  0.  Equation  (5)  is  finite  differenced  to 
obtain  the  sequence 


v'+a,/2  =  v'-a"2  +  A/ 


v  ^  m  J 
t  +  &irz 


r  =  r  +Atvr  ,  and 
9  r'+A/ 


(6) 


The  latter  guarantees  conservation  of  angular  momentum. 
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In  spherical  coordinates,  the  equation  of  motion  becomes 


dv,  m(ve2  +  v2) 
m  —  =  Fr  4 - , 

d< j> 

~dt 

d  (mr2  sin2  Qdtydt)  . 

- - =  r  sm  , 

dt 

where  r,  0,  and  $  are  the  spherical  coordinates  and  v#  =  r  sin  Qdtydt.  For  the  spherical  model, 
B  =  0  so  F0  =  F)  =  0  and  Fr  =  qEr.  If  dty/dt  =  0  initially,  then  it  remains  zero.  The  motion  is  then 
in  a  plane  <}>  =constant  through  the  polar  axis,  in  which  r  and  0  play  the  role  of  plane  polar 
coordinates.  The  angular  momentum  mrv8  is  then  constant,  and  Eq.  (7)  reduces  to  the  cylindrical 
case,  Eq.  (5). 


d(rv6) 
m— — 
dt 


=  rFa+mr2  sin  0  cos  0 


II.  POTENTIAL  EQUATION 


0  / 

\ - 1 — x 


Figure  3.  Configuration  in  planar  coordinates  with 
series  RLC  circuit  and  voltage/current  source.  In 
these  coordinates,  the  particles  are  charge  sheets 
with  motion  in  the  x-direction. 


The  planar,  cyclindrical,  and  spherical  configurations  for  one-dimensional  bounded  plasma 
systems  are  shown  in  Figures  3  and  4.  The  current  in  the  external  circuit  interacts  with  the  plasma 
current  via  surface  charge  on  the  electrodes.  Similarly,  the  potential  within  the  plasma  region  is 
affected  by  the  distribution  and  motion  of  space  charge,  the  electrode  surface  charge,  and  the 
current  in  the  external  circuit.  Thus,  we  seek  a  simultaneous  solution  for  the  potential  and  circuit 
equations. 
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Figure  4.  Configuration  in  cylindrical  and  spherical 
coordinates  with  series  RLC  circuit  and  volta¬ 
ge/current  source.  Particles  are  infinite  annuli  with 
motion  in  the  radial  direction  in  cylindrical 
coordinates,  and  spherical  shells  with  radial  motion 
in  spherical  coordinates. 


The  boundary  conditions  for  the  potential  equation  are  obtained  by  applying  Gauss’  Law 
to  the  system. 


r  Co  A+CT+  +  A_a_ 

jErfS=J^V'+-^I - =  0,  (8) 

S  V 

where  the  surface  S  encloses  the  plasma  and  electrodes.  A+  refers  to  the  surface  area  of  the 
positively  biased  (left/inner)  electrode,  A.  to  the  negatively  biased  (right/outer)  electrode,  and  a 
is  the  surface  charge  on  the  respective  electrode.  Note  that  p  has  units  of  charge/volume  and  a 
has  units  of  charge/area.  Equation  (8)  is  a  statement  of  Gauss’  Law;  the  first  part  reflects  the 
assumption  of  an  ideal  conductor  connecting  the  electrodes  to  the  external  circuit  elements,  the 
second  part  expresses  conservation  of  charge  in  the  system. 


Applying  Gauss’  Law  about  each  node  of  the  gridded  system,  and  using  the  definition  of 
potential,  we  obtain 


planar 


Ax2  e 


cylindrical 


r.Ar 

rj  +  ♦ 1  -  2 +  rj .  m<&j  _ ,  =  — ~  p j ,  and 


(9) 


spherical 


'>  +  \I7 


i^V+i  (rj  +  in  +  rj 


.  ..2  (r/+ 1/2  ~ r]- 1/2)^  _ 

-  \nfoj  +  rj  -  - 1  ~  P> 
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In  the  cylindrical  and  spherical  forms,  rJ±xa  =  rJ±Arll.  For  all  Eq.  (9),y=l,  2, ...,  nc- 1,  where  nc 

is  the  number  of  cells  in  the  gridded  space.  These  results  are  equivalent  to  the  flux-conserving 
method  of  Birdsall  and  Langdon  [2].  The  planar  and  cylindrical  result  of  Eq.  (9)  can  also  be 
obtained  by  applying  a  central  difference  to  the  Poisson  equation;  the  Finite  difference  result  is 
different  in  spherical  coordinates. 

For  a  one  dimensional  system,  the  boundary  conditions  can  be  written 

<1^=0  and  (10) 


„  (11) 

£.=7- 

Equation  (10)  fixes  a  reference  potential  for  the  system  without  implying  a  grounded  electrode. 
For  the  cylindrical  and  spherical  models  the  inner  electrode  is  driven;  the  outer  electrode  serves 
as  the  reference  potential  for  the  system  even  when  the  inner  electrode  is  not  present.  Equation 
(11)  can  be  written  at  one  half  grid  cell  from  the  boundary  in  conjunction  with  a  central  difference 
applied  to  the  definition  of  potential  to  obtain 


planar 


cylindrical 


spherical 


^1/2  ~ 


if  Ax 

^-%K+p»t 


<b0  i  [  r0  Pof  ro  1 1  j 


F  = 


if  rl  p0 


-i  °*A+3f»-A 


Equation  (9)  and  its  boundary  conditions  for  the  gridded  system  are  written  in  a  general 
matrix  form. 
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The  superscript  indicates  the  quantity  is  evaluated  at  time  t.  The  matrix  elements  in  planar 
coordinates  are 


b0  =  -l. 


Po 

Ax  +  2 


sP 

II 

<N 

II 

..,nc  -  1  ; 

b,  =  -2. 

ri 

II 

.  .,nc  -  1  ; 

II 

' — 

li 

© 

»— • 

..,nc  -  2  ; 

CL 

II 

■xT 

y = 1,2,. 

..,nc  -  1  ; 

.  Ax2 
1  ~  e 


The  matrix  elements  in  cylindrical  coordinates  are 


a;  ~  rj-V2' 

7  =  1,2,. 

,,nc  -  1  ; 

_  rv2  ’ 

bJ  =  ~2rj  > 

7  =  1,2,. 

.  ,,nc  —  1  ; 

Cj  =  rj  +  1/2  > 

y  =  0.1,. 

.,nc  -2  ; 

,  r 0  Po  2  2, 

~&+Ar  2Ari>m  o)  ’ 

4  =  rjPj  > 

7-1,2,. 

..,r,c  - 1  ; 

The  matnx  elements  in  spherical  coordinates  are 


K  —  r\n  > 


ai  ~  rj  -  1/2  > 

7  =  1,2,. 

..,nc  -  1  ; 

by  =  ~ {^y  -  1/2  rj  +  1/2)  > 

7  =  1,2,. 

..,/ic  - 1  ; 

Cy  =  ry  +  1/2  ’ 

7=0,1,. 

. . .  nc  —  2  : 

(16) 

^y  =  (ry  +  1/2  _  ^y  -  l^)Py  > 

7  =  1,2,. 

..,/ic  - 1  ; 

A r 

and  3i- 

When  the  center  conductor  is  not  present  in  the  curvilinear  models,  the  boundary  conditions 
must  be  modified.  From  Gauss’  Law,  the  electric  field  at  the  origin  must  be  zero.  Integrating 
Gauss’  Law  from  the  origin  to  r  -  rxa,  we  obtain  the  modified  form  of  Eq.  (12)  for  the  hollow 
cylindrical  system  (2,  Section  14-10], 
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CIRCUIT 

The  external  circuit  is  coupled  to  Eqs.  ( 1 3)-(  1 6)  through  conservation  of  charge  at  each  wall, 

AAa-  Qcoxv  +  AQ  ,  (20) 

where  Qcoin  is  the  charge  deposited  by  the  convection  (particle)  current  and  A Q  is  the  charge 
deposited  by  the  external  circuit  current,  both  over  some  interval  in  time.  Equation  (20)  is  applied 
at  the  positively  biased  electrode  as  shown  in  Figures  3  and  4,  guaranteeing  conservation  of  charge 
at  all  times.  The  same  logic  can  also  be  applied  to  the  other  electrode;  however,  the  surface  charge 
on  the  second  electrode  is  determined  readily  from  Eq.  (8)  when  the  first  surface  charge  is  known. 
The  charge  conservation  equation  becomes 


o'  =  o'A'  + 


QL^  +  Q'-Q' 


where  Q  is  the  charge  on  one  plate  of  the  external  circuit  capacitor.  An  alternate  method  of 
coupling  the  circuit  to  the  potential  matrix  is  applying  continuity  of  current  (Kirchhoff  s  Current 
Law)  at  the  boundary  [2,  Section  16-9], 

to  ,L  (22) 

dt  A  ’ 

where  Jcomv  is  the  plasma  convection  current  at  the  electrode.  The  methods  are  equivalent  when 
a  first-order  backward  difference  is  used  for  dc/dt  and  /  =  dQldt.  Since  Qconv  is  in  general  a  noisy 
quantity  in  a  particle  simulation,  any  other  quantity  in  Eqs.  (20)  and  (22)  will  contain  similar  noise. 
Thus  Eq.  (20)  causes  the  wall  charge  a  to  be  noisy  as  might  be  expected,  because  the  capacitor 
charge  reacts  to  the  particle  convection  current  only  through  the  wall  charge;  i.e.,  particles  absorbed 


8 


by  the  wall  contribute  immediately  to  c,  but  the  charge  drains  slowly  to  the  capacitor  through 
currents.  It  can  be  shown  that  Eq.  (22)  results  in  the  convection  current  being  absorbed  gradually 
into  a,  so  the  noise  is  induced  in  the  capacitor  charge  Q  (and  consequently  in  the  external  current 
[)  to  satisfy  conservation  of  charge.  Therefore  we  use  the  conservation  of  charge  method  of  Eq. 
(20). 


III.  GENERAL  SERIES  RLC  CIRCUIT 

Four  cases  cover  the  full  range  of  external  circuit  parameters.  For  the  general  voltage-driven 
series  RLC  circuit,  the  capacitor  charge  Q  is  advanced  using  Kirchhoff’s  Voltage  Law, 


L 


SQ 

dt 2 


+  R^~  +  ®=v(t)  +  <t>«-4>o 


(23) 


The  polarity  of  the  source  and  resultant  positive  current  are  shown  in  Figures  3  and  4.  The 
general  circuit  equation  is  finite  differenced  using  the  second-order  backward  Euler  representation 
of  the  first  derivative, 

'dQ  Y_  3Q' -4Q'-*  +  Q'~2a‘  (24) 

K~dt  J  ”  2Ar  ’ 

and  the  second  derivative, 


VgY  3(dQ/dt)‘ -4(dQ/dt)"“  +  {dQ/dt)‘-2A‘ 
y  dt 2  j  2 Af 


9Q'  -  24Q'~M  +  22Q'~2A;  -  8Q'~,a;  + 
4A t2 


(25) 


The  latter  is  obtained  by  a  second  application  of  the  first  derivative  to  Q.  An  alternate  4  point 
difference  for  the  second  derivative  is  given  by: 


V<2  Y  _2Q‘  -5Q"‘y  +  4Q"^ -Q'-^ 
y  dt 2  J  A/2 


(26) 


The  charge  on  the  capacitor  is  not  known  at  t.  Combining  Eqs.  (23)-(25),  we  obtain 


Q‘  = 


V(t)  +  <V'C-<V,-K' 

Oo 


(27) 
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where 


K‘  =  a,Q"*  +  0,0' +  0,0' - +  a<Q 


I  -4A( 


9  L  3  R  1 
a°~4^  +  2  *  +  C’ 
.  L  R 

ai  =  -^7i-2  — , 

A/  A/ 

11  L  1  /? 

2  A/:  +  2A/’ 

(X3  =  -2  — r ,  and 
Ar 
1  t 

0C4  =  - — r  . 

4  a/2 


(28) 


Combining  the  potential  equation  results,  Eqs.  (13)-(16),  with  the  circuit  equation  results, 
Eqs.  (27)  and  (28),  using  the  boundary  condition,  Eq.  (21),  we  obtain  the  self  consistent  field 
solution  matrix  for  the  voltage-driven  series  RLC  circuit  case.  The  matrix  can  still  be  represented 
by  Eq.  (13),  replacing  elements  of  Eqs.  (14)-(  16)  as  follows.  In  the  planar  model, 


bn  =  - 1  and 


,  Po  Of*  1 

dn  =  —  +  —  +  ■ 


0  2  Ax  -4  Ax 

In  the  cylindrical  model, 


0 idEA 

V  «<> 


Ar 


b°=-r'°-^  and 


d„  = 


2  2 
rin-r«  .  r{ 


1/2  '0  .  '0  ,-etj  1 


r 


0  2Ar  Po  +  Ar°*  +  2nhAr^conv  ®  ? 


.  V(0-AT 


V 


«o 


In  the  spherical  model, 


(29) 


(30) 


,  2  A r 

bo  =  -r\a~- -  and 


4JKCL0 


4>  =  (ria  ~  ro  )Po  +  3r0  ^ 


QL.-Q'~*  + 


1  - ai  .  V(r)-/C 


Oo 


(31) 
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Here,  A  is  area  of  the  planar  electrodes  and  h  is  axial  length  of  the  cylindrical  system.  The  solution 
is  then  self-consistent  and  second-order  accurate  for  the  general  circuit  case.  The  matrix  can  be 
solved  using  any  algorithm  optimized  for  tridiagonal  matrices  [3]. 

IV.  OPEN  CIRCUIT  (FLOATING  OUTER  ELECTRODE) 

When  C  -» 0,  the  impedance  of  the  external  circuit  approaches  infinity,  becoming  an  open 
circuit.  The  potentials  on  the  boundaries  are  floating;  no  circuit  solution  is  required  since  there 
is  no  external  current.  The  surface  charges  on  the  electrodes  influence  the  potential  as  always, 
but  the  electrodes  cannot  exchange  charge  via  external  current.  In  this  case,  the  field  solution  is 
given  by  Eqs.  (13)-(16),  with 


(32) 


V.  SHORT-CIRCUIT 

When  R-L=0  and  C  -» the  external  circuit  is  a  short,  with 

The  short-circuit  case  is  applied  in  practice  when 


(33) 


planar 

cylindrical 

spherical 


C  >  10s , 


eA/l 

C 


2n£hl\n(rK/r0) 

C 


4ner«ro /(r*-r0) 


>105, 


>  10s, 


(34) 


where  /  is  the  length  of  the  planar  plasma  region. 

The  field  solution  is  still  given  by  Eqs.  (13)-(16),  with 

a\  =  ^0  =  ~  ^0  =  ®  » 
and 


(35) 
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planar 


cylindrical 

spherical 


,  .n,  V(Q 
Pi  j  > 

J  ,  riaV(t) 
d\  ~  riPi - j — .  and 

.  3  >. ,  i&m 

d\  ~  (r3 a  r\n)V\  y 


(36) 


Equation  (35)  eliminates  the  first  row  of  Eq.  (13).  In  Eq.  (36), /depends  on  the  model  as 
given  in  Eqs.  u4)-(16).  Note  that  the  wall  charge  is  no  longer  required  to  solve  the  potential 
equation.  Wall  charge  is  found  using  Eq.  (12)  once  the  potentials  have  been  determined,  and  the 
current  is  found  by  finite  differencing  Eq.  (22), 

&-&-*)  (37) 

Ar  )‘ 


Determining  the  current  in  this  way  produces  a  noisy  result  as  discussed  above;  however,  with  a 
short  between  electrodes,  we  expect  large  currents  with  rapid  changes  since  potential  differences 
cannot  exist  along  an  ideal  conductor.  Note  that  here  /  is  only  a  diagnostic  quantity,  so  the 
time-centering  is  not  a  problem. 


VI.  CURRENT-DRIVEN  CIRCUIT 

The  final  case  is  the  current-driven  external  circuit.  An  ideal  current  source  is  assumed 
which  can  drive  the  specified  time-varying  current  /(f).  The  external  circuit  elements  R,  L,  and 
C  are  ignored  since  an  ideal  current  source  is  an  open  circuit.  Then  Eqs.  (13)-(16)  are  applied 
with  the  wall  charge  found  by  finite  differencing  Eq.  (22)  for  diagnostic  purposes. 


VII.  INITIAL  CONDITIONS 

The  multi-point  finite  difference  methods  require  initial  values  for  the  Q ",  where  n  <  0. 
Physically,  these  values  are  used  to  obtain  the  desired  initial  conditions  for  the  circuit  equation, 
Eq.  (23).  For  example,  the  initial  charge  on  the  capacitor,  Q°,  and  the  initial  current  in  the  external 
circuit,  f,  form  a  complete  set  of  initial  conditions  for  the  differential  equation.  However,  the 
finite  difference  requires  five  initial  values  for  Q  (four  for  the  4-point  method).  There  are  several 
ways  the  conditions  can  be  obtained. 

The  traditional  method  for  starting  a  multi-point  scheme  (second  or  higher  order  accurate) 
is  to  use  a  2-point  method  (first  order  accurate)  to  obtain  the  required  initial  values.  A  smaller 
timestep  is  used  with  the  2-point  method  to  maintain  the  same  accuracy.  This  presents  a  problem 
for  a  PIC  code;  the  time-centered  mover  is  initialized  such  that  positions  are  known  at  integral 
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timesteps,  while  velocities  are  known  at  half  timesteps  [2].  Thus,  it  is  difficult  to  switch  to  a  new 
A t  and  maintain  the  time  centering.  Also,  switching  schemes  is  inefficient  from  a  coding  stand¬ 
point.  In  addition,  the  stability  of  the  starter  method  must  be  considered  in  relation  to  the  circuit 
parameters  R,  L,  and  C. 

Another  method  of  initializing  the  solver  is  to  solve  the  circuit  equations  analytically.  To 
do  this,  we  must  replace  the  plasma  by  a  known  impedance.  Using  the  vacuum  capacitance  of 
the  plasma  region  is  the  obvious  choice;  physically,  this  means  there  is  no  plasma  until  t- 0+.  If 
plasma  is  then  introduced,  the  impedance  changes  abruptly  and  the  circuit  has  been  conditioned 
for  a  different  system.  This  problem  is  less  severe  when  the  plasma  is  generated  at  a  slow  rate 
since  the  impedance  change  is  gradual. 

If  the  method  turns  out  to  be  stable,  the  initial  conditions  will  be  damped  regardless  of  the 
value  (this  includes  desired  initial  conditions  as  well  as  error  in  the  initial  conditions).  If  the 
method  is  unstable,  any  error  in  the  initial  conditions  grows  exponentially.  If  the  method  is 
marginally  stable,  any  error  in  the  initial  conditions  remains  in  the  solution,  neither  growing  nor 
damping. 

VIII.  STABILITY 

We  now  explore  stability  of  the  circuit  equation,  Eq.  (27).  As  is  customary  for  stability 
analysis  [4],  we  neglect  the  driving  terms  and  study  the  homogeneous  circuit  equation 


dt2  dt  C 


(38) 


We  study  the  stability  of  the  5-point  circuit  difference,  Eqs.  (24)-(25),  as  well  as  the  4-point 
difference,  Eq.  (24)  and  (26). 


In  the  limit  of  no  inductance,  L  0,  both  methods  produce 


Q‘ 


2Ax 


i-M  .  ,->1-2^ 


3+«c  re  +Q 


=  o 


(39) 


Letting  Q‘  =  Q°ev  and  4  =  e1’4'  we  obtain 

Q'  =  4Q'"A'  =  42C?"M;  -  (40) 

where  |  4 1^  1  is  required  for  stability.  Here,  y  and  4  316  arbitrary  complex  variables.  Then  the 
characteristic  stability  equation  for  Eq.  (39)  is 

42(3  +  2A///?C)-44+l=0.  (41) 
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The  roots  are 


2  +  V1-2A t/RC  (42) 

3  +  2A tIRC  ' 


Rgure  5.  Stability  roots  in  the  limit  L  ->0.  Since 
|  %  £  l  everywhere,  the  method  is  stable.  The 
scheme  can  only  follow  the  RC  time  when  to  <  RC/2. 


As  shown  in  Figure  5,  both  methods  are  stable  in  the  limit  L  -»  0  for  all  positive,  real  At/RC. 


Now  we  attack  the  more  difficult  general  case.  The  general  characteristic  stability  equations 
for  the  four  and  five  point  schemes  are  respectively 


^2  +  |t1  +  ^-42(5  +  2t,)  +  4  4+|x, 


-1=0  and 


£*(9  +  6x,  +  4x^)  -  43(24  +  8x.)  +  ^2(22  +  2t,)  -  +  1  =  0  , 


(43) 

(44) 


where  the  normalized  times  are  x,  =  R  Ail L  and  Xj  =  AtHLC.  We  obtain  the  roots  of  Eq.  31  using 

the  Lin-Bairstow  method  [5],  which  gives  the  complex  roots  of  polynomials.  Figures  6  and  7 
show  the  stability  of  the  four  and  five  point  methods,  respectively,  for  a  wide  range  of  x,  and  x2. 
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(C)  (d) 

Figure  7.  Magnitude  of  the  four  stability  roots  of  the  five  point  method,  whose  characteristic 
equation  is  Eq.  (44).  Since  all  four  roots,  |  £1iJ>3i4  |<  l,  the  five  point  difference  scheme  is  stable 

over  the  range  shown. 


IX.  SIMULATION  ACCURACY 

As  discussed  above,  the  circuit  and  Field  solution  is  second-order  accurate.  We  now  dem¬ 
onstrate  the  accuracy  of  the  five  point  circuit  as  implemented  in  the  code,  PDP1 .  To  compare  the 
simulation  results  with  analytic  circuit  theory,  the  permittivity  of  the  plasma  region  is  taken  as 
1020  and  no  plasma  is  used.  Then  we  have  a  passive  voltage-driven  series  RLC  circuit.  The  current 
for  a  sinusoidal  driving  voltage  V  =  sin(ow  +  9)  can  be  predicted  using 


a2V/(coZ)  cos(0  -  5)  +  VIZ  sin(0  -  5) 

' " - ssrn - exp(a',) 


a,W(o)Z)cos(0-6)  +  V7Zsin(0-8)  V 

+ - - - t - exp(a,r ) + — sin(cor  +  0  -  5) , 

flj/ d2  1  Z 


where  co  is  the  angular  frequency  of  the  voltage,  0  is  the  initial  phase  of  the  voltage,  and 
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Figure  8.  The  frequency  spectrum  of  the  steady-state  current  obtained  from  PDP1 . 

We  choose  R  =  1,  C  =  5  x  10-6,  L  =  KT6,  and  0)=  106,  and  initial  conditions  Q(0)=0  and  /(0)=0. 
The  code  PDP1,  using  the  same  parameters  and  timestep  At0  =  27t/128co,  gives  the  results  shown 
in  Figures  8  and  9.  Here,  T,  =  0.049  and  r2  =  0.022  for  the  baseline  case.  Note  that  the  t  scale 
with  Ar/Af0. 

From  the  frequency  spectrum  of  the  current  shown  in  Figure  8,  we  see  the  driven  frequency 
peak  is  six  orders  of  magnitude  larger  than  the  magnitude  at  other  frequencies,  indicating  nearly 
a  pure  sinusoid.  Since  PDP1  stores  results  in  single  precision  (32  bits),  we  expect  roundoff  error 
in  the  sixth  or  seventh  significant  digit,  so  the  powers  less  than  10'7  are  neglible.  Comparing  the 
phases  I(t),  we  see  the  the  PDP1  results  follow  Eq.  (45)  closely,  with  increasing  phase  error  as  At 
increases.  From  the  historv  of  the  current,  we  see  the  initial  transient  due  to  the  charging  of  the 
capacitor  from  Q( 0)=0. 
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The  relative  error,  plotted  in  Figure  10,  follows  the  curve  0.001 58(A/)' 94  closely.  An  ideal 
second  order  accurate  scheme  would  result  in  a  power  fit  exponent  of  2.  This  demonstrates  second 
order  accuracy,  with  errors  resulting  from  truncation  of  the  finite  difference  at  At 2  terms. 

X.  CONCLUSION 

A  method  for  the  simultaneous  solution  of  the  coupled  potential  and  external  circuit 
equations  for  one-dimensional  electrostatic  plasma  particle  simulations  is  presented.  The  method 
is  stable  over  many  orders  of  magnitude  for  the  values  of  the  RLC  circuit  elements,  and  can  in 
principle  be  extended  to  arbitrary  external  circuits. 

The  method  is  implemented  in  the  codes  PDP1  (Plasma  Device  Planar  1  Dimension),  PDC1 
(Cylindrical),  and  PDS1  (Spherical)*.  These  codes  have  been  used  to  simulate  many  complete 
bounded  plasma  devices  [6-11],  including  voltage-driven  RF  discharges,  plasma  immersion  ion 
implantation  devices,  and  Q-machines.  The  codes  have  performed  reliably,  generating  many 
interesting  discussions  and  discoveries. 
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